clear; clc;
% 在不同 E (0.2, 1.2, 5) 下
% 膜电位 x(t) 和能量 H(t) 的时间序列，并在标题中显示后期平均能量 <H>

% 参数设定
a     = 0.7; 
b     = 0.8; 
c     = 0.1; 
alpha = 0.25; 
beta  = 0.01; 
delta = 0.1;

% 外部激励幅值和频率
A_val = 0.35;
f_val = 0.72;

y0 = [0.2; 0; 0.01]; % 初始条件

% 积分时间设定
T_total = 1200;  % 总模拟时间
T_trans = 400;   % 舍去瞬态时间，后 800 时间单位用于计算 <H>
h       = 0.001;  % RK4 步长

% 待模拟的 E 列表
E_list = [0.2, 1.2, 5];

midnightPurple = [0.2, 0, 0.4];
cyanGreen      = [0, 0.8, 0.6];

tlo = tiledlayout(3,2);

for i = 1:3
    E_val = E_list(i);

    % 将参数打包，含 E 字段
    params = struct('a', a, 'b', b, 'c', c, 'alpha', alpha, ...
                    'beta', beta, 'delta', delta, ...
                    'A', A_val, 'f', f_val, 'E', E_val);

    % a) 数值积分 (使用 hybridNeuron_E)
    [t, Y] = rk4(@(tt,yy) hybridNeuron_E(tt, yy, params), [0, T_total], y0, h);

    % 提取 x, y, z
    x_data = Y(:,1);
    y_data = Y(:,2);
    z_data = Y(:,3);

    % b) 根据论文公式 (7) 计算能量 H(t)
    %   H = 0.5*x^2 + (1/(2c))*y^2 + 0.5*b*z*sin(beta*z)*y
    H_data = 0.5 * x_data.^2 ...
           + (1/(2*c)) * (y_data.^2) ...
           + 0.5 * b .* z_data .* sin(beta .* z_data) .* y_data;

    % c) 除瞬态 (t > T_trans)，计算后期平均能量
    idx = (t > T_trans);
    t_plot = t(idx);
    x_plot = x_data(idx);
    H_plot = H_data(idx);

    mean_H = mean(H_plot);  % 后期平均能量

    % d) 绘制每个 E 下的 x(t) 与 H(t) (3行2列)

    % 左侧子图：x(t)
    nexttile(2*(i-1)+1);
    plot(t_plot, x_plot, 'Color', midnightPurple, 'LineWidth', 1.5);
    xlabel('\tau'); ylabel('x');
    title(sprintf('E=%.1f', E_val));
    grid on;

    % 右侧子图：H(t)
    nexttile(2*i);
    plot(t_plot, H_plot, 'Color', cyanGreen, 'LineWidth', 1.5);
    xlabel('\tau'); ylabel('H');
    title(sprintf('E=%.1f, <H>=%.5f', E_val, mean_H));
    grid on;
end